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1. Introduction 

The numerical evolution of Einstein field equations has proved to be a formidable 
task. Ill-posed formulations of the dynamical equations, the presence of singularities, 
exponential growth of constraint violations, inadequate inner (in the case of black 
hole excision) and outer boundary conditions, and the lack of robust shock-handling 
hydrodynamical algorithms for strong gravity regimes are some of the open problems 
encountered in numerical relativity. To this list, we should add the problems intrinsic 
to the generation of astrophysically realistic initial data (ID) sets, such as the lack of 
adequate formulations to describe accurately the gravitational radiation content at the 
initial time step, the difficulty in specifying the fluid motion corresponding arbitrary 
stellar spins, or the inner boundary conditions around black hole singularities, to 
mention just a few. While none of the above issues has been solved to satisfaction, 
great progress has been made in the past decades in all these directions. 

Controlling the exponential growth of constraint violating modes is essential for 
the stability of numerical simulations. Some of these modes are produced in the bulk of 
the numerical grid, originating in small violations due to roundoff and truncation errors 
which are inevitable in any numerical treatment. Other modes are generated at the 
grid boundaries by inappropriate boundary conditions (see, for instance, [I1I21)- While 
there is a purely numerical component to these instabilities, the scientific community 
consensus is that the choice of dynamical formulation (a given way of casting the Einstein 
field equations) will play an important role in their control |S1 IH E]- Most of the 
earliest attempts at modeling numerically general relativistic systems were based on 
a formulation introduced by Arnowitt, Deser, and Misner j^j, widely known as ADM. 
However, it soon became clear that the life of the simulations could be extended by 
modifying ADM or simply replacing it with different formalisms. A method originally 
developed by Shibata and Nakamura [Zj and later used by Baumgarte and Shapiro 
[Hj (BSSN) is today the most commonly used for three-dimensional simulations. The 
literature offers an ever-growing list of new evolution formulations which can be divided 
in two groups: unconstrained [3011111113110111111121113111111311311111111113 and 
constrained [23I2II- Some of these have been tested numerically under conditions that 
are either easy to implement numerically and / or have a high degree of symmetry. 
Test cases have been constructed based on perturbations of Minkowski spacetime, 
linearized waves, non-linear plane waves, Gowdy and Robert son- Walker metrics. Brill 
and Klein-Gordon waves, non-linear versions of Maxwell equations, and single black hole 
spacetimes (for a comprehensive analysis of such tests, we refer the reader to j22j). The 
importance of these tests is twofold: they show the weaknesses and strengths of each 
formulation in easy to understand cases and they serve as a first discriminatory round of 
benchmarks that are simple enough to be in the reach of most research groups. However, 
binary simulations have been performed sparsely: they require very complex numerical 
codes as well as non-trivial computational power. To our knowledge, ADM and BSSN 
are the only formulations that have been tested in three-dimensional finite-size compact- 
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object binary simulations: either with neutron stars [231 1211 12^1 12ni I2Z1 12H1 121] or black 

holes [3niEIlE21Ell311- 

In this paper, we design and implement a numerical setup based on binary 
neutron stars (BNS) simulations. BNS inspiral scenarios are fully general relativistic, 
three-dimensional, and do not present singularities, thus decoupling the problem of 
inner boundary conditions needed for black hole excision. They provide a robust 
test for the capacity of the numerical implementation to handle angular momentum 
conservation, usually one of the most problematic quality control monitors f. They also 
are astrophysically realistic systems which are related to some of the most important 
observational phenomena of our times: gamma-ray burst engines and generators of 
gravitational waves. BNS simulations usually require extensive computational resources 
and the length of the runs could, in principle, render these tests impractical. Here we 
show how small, low resolution grids can be used to gain insight into the stability of 
different numerical schemes, with runs that only take a few hours on single-processor 
workstations. 

To illustrate the efficacy of the proposed testing ground, we introduce a new 
technique that was developed with the help of short trial-and-error BNS simulations. 
This method is based on the approximate satisfaction of the Hamiltonian constraint, 
obtained by relaxing the conformal factor. Anderson and Matzner [21] recently 
performed black-hole simulations where the Hamiltonian and momentum constraints 
elliptic equations are solved at every time step. In our approach (which for this 
paper involves only the Hamiltonian constraint) the conformal factor is relaxed from 
a parabolic equation that drives the Hamiltonian constraint into exponentially decaying 
solutions. This method is an adaptation of the K-Driver algorithm introduced by 
Balakrishna et al. [33], used to enforce the maximal slicing condition. This technique 
reduces the Hamiltonian constraint violation by a factor of five with respect to the ID 
values, and by two orders of magnitude with respect to BSSN simulations. Berger j3£] 
has shown recently how the satisfaction of the Hamiltonian constraint throughout the 
evolution plays a role in the overall quality of the simulation. We observe a similar effect 
in our runs which is particularly evident in the remarkable improvement of the behavior 
of the total angular momentum of the binary. The total angular momentum as a function 
of time agrees with post-Newtonian estimations for more than one orbital period for 
the case of BNS outside the innermost stable circular orbit (ISCO). In addition, the 
computational overhead caused by the relaxation method adds only an extra 5% to the 
running time of the simulations, making the technique more practical than solving the 
elliptic equation derived from the Hamiltonian constraint. 

Section |21 describes the time evolution equations for the gravitational fields, 
summarizing the differences between BSSN and the Hamiltonian relaxation method 
with special emphasis in the boundary conditions employed in this paper. Section |H1 

f Total angular momentum is not strictly conserved in general relativity. However, most of the 
variations observed in current simulations are due to numerical error and not to the physical loss 
due to gravitational radiation. 
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describes the BNS numerical experiment, introducing the new full general relativistic 
hydrodynamical code (GRHyd) employed for the time evolution. Section 0] presents 
results obtained with the Hamiltonian relaxation method and compares them with 
the corresponding BSSN counterparts. The Appendix section shows validation and 
convergence tests, as well as the derivation of the post-Newtonian (PN) calculation for 
the loss of angular momentum for point-mass binaries. 



2. Equations for the Gravitational and Hydrodynamical Fields 

2.1. Time Evolution Equations 

In this paper we use geometrized units {G = c = 1) and the Greek (Latin) indices run 
from to 3 (1 to 3). 

Numerical relativistic treatments usually cast the metric in the "3-1-1" form 

ds^ = -a^dt^ + j,j{dx' + (3'd t){dx^ + (3^dt) , (1) 

where a, P\ and 'jij are the lapse function, shift vector, and spatial metric tensor, 
respectively. The extrinsic curvature Kij is defined as 

Kij = -{dt - CpH, I (2a) , (2) 

where £^ is the Lie derivative with respect to /3* 0. Using these fields, we can rewrite 
Einstein's field equations 

G^y = SttT^,^ (3) 
set of four differential equations: 

R - KijK'^ + = levrp , (4) 
DjKi - DiK = SnSi , (5) 

known as the Hamiltonian and momentum constraints, and 

{dt - = -2aKij , 

{dt - Cp)K,j = -DiDja + a{Rij - 2KuK'j + KK.j - %Ti[Sij + ^7,,(p - S)]} , (6) 

which provide the evolution in time of the spatial metric and the extrinsic curvature. 
The symbol Di represents the covariant gradient with respect to the tensor 7jj. The 
fields p, S", and Sij are derived from the matter fields by sphtting the stress-energy tensor 
T^u in components parallel and perpendicular to the normal of the spatial hypersurface 

Si = -lic.npT'^^ , 

Sij lialj^T , 

S = I'^S,, . 

These equations are the basis of the ADM formulation [H] . 
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Following York 3?^, we rewrite the metric and the extrinsic curvature as 

These decompositions define the fields if), '^ij, Aij, and K, known as the conformal factor, 
the conformal metric, the conformal traceless extrinsic curvature, and the trace of the 
extrinsic curvature respectively. Using these variables, the Hamiltonian and momentum 
constraints can be rewritten as 

fW,D,^p -^R+ ^A,A'^ - ^K' + 2n^'p = 0, (7) 

Djitp'^A^') - -tp^D'K - S-K-^^S' = . (8) 
3 

The ID sets are generated by solving (usually numerically) these equations. These 
solutions are correct within the bounds of the truncation error associated with the 
finite-difference scheme of choice. These initial constraint violations will be amplified 
when using unconstrained formulations like ADM and BSSN. 

Equations for the time evolution of the fields ^ij, K, and Aij can be derived from 

®: 

{dt - = -2aAij 

{dt - Cp)K = -fWjDia + -aK'^ + aAi^A'^ + A-na{p + S) 

3 

{dt - Cp)A,j = ij-\-D,D,a + a{R,j - SvrS,,)]^^ + a{KAij - 2AaA' ,) , (9) 

where the superscript TF indicates the trace-free part of the tensor. These fields are 
complemented with the variable known as the conformal connection, introduced in [7j 

r = -f ^„ (10) 

where we follow the notation of [8^ . An evolution equation is derived for these variables 
from © and ©: 

dtV = dj{2aA'^ + Cpf^) 

- 2a {h'^'K,, - QA'^ [ln{tp)ij - r.^A^^ + Mf^S^) . (11) 
Equations Q and (jllj) together with a time evolution equation for the conformal factor 
{dt-Cp)\n{,l,) = -\aK (12) 



are the basis of the BSSN formulation [HIHI- 

We define the Hamiltonian constraint residual Ti as the r.h.s. of equation (|7j). Note 
that Ti will only be null when dealing with exact solutions in the continuum case; in 
any numerical treatment the round-off and truncation errors will violate the constraint 
((Zj). This equation is an elliptic second-order PDE that, in principle, can be used to 
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determine the conformal factor. Anderson and Matzner |2II performed simulations of 
static single black-hole spacetimes, solving equation (|7j) at every time step. 

The Hamiltonian relaxation technique proposed here derives the conformal factor 
from solving approximately the Hamiltonian constraint, in a way that resembles the K- 
Driver algorithm [HSI used for calculating the lapse function. Instead of solving equation 
([7j), we compose an alternative parabolic equation which will relax to a solution of 
the Hamiltonian constraint through an iterative scheme. An example of such equation 
is 

dt'i^ = en {deU + Vh H) , (13) 

where en and rjn are fine-tuning parameters. Note that the time coordinate t' is not 
the physical time, but a relaxation parameter that has only numerical meaning: for 
every At of physical time, equation (fT^ will be relaxed up to a maximum number of 
times At /At' (see below). A stationary solution (in t') of equation ()13|1 enforces the 
condition dt'Ti. = —rjH 'H-, which corresponds to an exponential decay of the Hamiltonian 
constraint [HK] . 

The dt'li. term of equation p3p could in principle be derived from the time evolution 
equations Q. However, the resulting equation is quite long and difficult to handle 
numerically. A much simpler alternative is the use of the Forward- Time Centered- Space 
(FTCS) finite-difference expression 

where Ti represents the constraint violation at the relaxation iteration step m. 

The relaxation of is performed at each one of the steps of the ICN method (one 
Predictor and two Corrector steps) and is done after the corresponding update of the 
rest of the gravitational fields and before the update of the hydrodynamical fields (i.e.; 
p). The relaxation in the Predictor stage follows these steps: 

1) Initial update of ip: 

^^^ = ^n-i + At' enVHUn-i, (15) 

where the subscript n follows the physical time step. Field values from the previous time 
step (i.e.; ipn-i, 'Hn-i, etc.) do not carry any upper index since they are not "relaxed". 
This update is followed by a re-evaluation of the boundary conditions for ip^ which are 
explained in detail in section IT^ 



2) Subsequent updates of ip: Loop on index m from 1 to M 

C = C + Ai' " " + VH nr~' ) , (16) 



where the Hamiltonian residual 7i„™ ^ is evaluated as 

m—1 u m—l(„i,m—l ~ 



7^ = K"^-^(C"\%- (n-1), A, (n-i),ete.) . (17) 
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These updates of t/'™ are also followed by a re-evaluation of the boundary conditions, 
after which the L2 norm ||7i„'"||2 is calculated . The relaxation iteration is stopped 
when ||7i„™'||2 < H'^^n-ilb or when the maximum number of iterations steps M (chosen 
here to be 25) is reached. 

The same steps are followed in the first Corrector (second Corrector) stage, but 
replacing the field values from the previous time step with the corresponding updates 
generated in the Predictor (first Corrector) step. 

2.2. Lapse and Shift Equations 

The lapse function a and shift vector /3* (i.e.; the gauge fields) are related to the 
coordinate degrees of freedom of the theory of general relativity. The choice of lapse 
function controls the way the spacetime continuum is split into a foliation of spatial 
hypersurfaces, while the shift vector indicates how the spatial coordinates change from 
one hypersurface to the next. The gauge fields are usually chosen to gain stability 
during the simulation, thus, they depend on the choice of time evolution methods 
and the particular physical system under study. Previous BNS studies |2H1 have 
shown that the maximal slicing condition (namely, K = 0) implemented with the K- 
Driver algorithm jHSI is a robust choice for the lapse function when using the BSSN 
formulation. Note that for our simulations maximal slicing conditions of any form have 
a clear advantage over other alternatives, given that the ID sets are based on this 
condition for the lapse. This condition is based on the parabolic equation for the lapse 
function a 

dt'a = -e^ {dt'K + ri^K) , (18) 

where ea and r/a are free parameters. 

The same studies also showed that the F-Driver jHE] condition for the shift 
vector makes possible long term BNS evolutions. A similar parabolic relaxation method 
is used here, this time to relax the shift vector in such a way that it drives the BSSN 
variable F* exponentially to zero. The corresponding parabolic equation is 

dt>(3' = ep (9i,r + r/^r) , (19) 

where e^? and rjp are the corresponding free parameters. We also used in our simulations 
the condition known as (3 freezing, that leaves the shift vector unchanged (i.e., identical 
to its initial value) all throughout the evolution. 

Finally, all the simulations presented in this paper were performed in the frame that 
rotates with the binary (corotating), which has been proved to enhance the stability of 
the evolutions not only in generally relativistic EH] , but also in Newtonian binary 
simulations jSHl- 

2.3. Boundary Conditions 

The choice of outer boundary conditions can be crucial for the stability of any numerical 
scheme. We adopt Sommerfeld (radiative) boundary conditions for the conformal metric 
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7jj, the traceless part of the extrinsic curvature Aij, and the trace of the extrinsic 
curvature K, and we set = at the boundaries (see Duez et al. j^Sj for details). 
These conditions are used for both BSSN and Hamiltonian Relaxation runs. 



2.3.1. Boundary Conditions for the Conformal Factor ip One important difference 
between the BSSN and the Hamiltonian relaxation runs is the boundary condition 
adopted for the conformal factor tp. The BSSN runs use the same (Sommerfeld) 
boundary condition employed for the other gravitational fields. 

The Hamiltonian relaxation runs use a boundary condition that enforces the 
satisfaction of 7-^ = in its finite-difference form. Equation ((7j) can be expressed as 

fW.D.tfj = d,d,^ - f'^dk^ = , (20) 

where collects all the terms that do not depend on derivatives of tp. Let's consider 
first the case of the grid points in the middle of the cube faces (i.e.; excluding the 
edges of the grid). If our numerical grid is a Cartesian cube with A^^ points, the points 
corresponding to the surface x = Xmax are represented by the indices (A^, j, k) \. As an 
example of our boundary condition for the conformal factor, we will derive an expression 
for the boundary values ip{xmax) = i^N,j,k, based on the finite-difference approximation of 
equation at the point next to the boundary (i.e.; the point with indices (A^ — 1, j, k)). 
There, the differential operators involving partial derivatives along the x direction are 
expressed as 

'^N,j,k ~ ^N-2,j,k 



dxdyil)' 
etc, 



2/\x 

n '^N,j,k ~ '^i'N-l,j,k + V'jV-2,j,A: 

(Axy 

'^N,j+l,k ~ '^N-2,j+l,k ~ '^N,j-l,k + V'jV-2,j-l,fc 

4 Ax Ay 



where n is the time step index and Ax and Ay represent the grid spacing along the x 
and y axis, respectively. Replacing the differential operators of equation (plj) with their 
finite-differencing counterparts results in an algebraic equation from which we obtain 
an expression for ip"^ ■ ^ of the form 

- i^iv.>(r) (Ax)^ 

\lxx)N-l,j,k - J- N-l,j,k 

where -FArj,A;(V'") represents an algebraic function of V'l^j,^; V^Tvj-i.fc; i^Nj+i^ki ''pN,j,k-v 
and ip%j^f,_^_i. Equation EH) defines a system of linear equations on the unknowns ipN,j,kJ 
whose solution can be very time consuming. A faster alternative is to evaluate the 
equation ^ using the values ^j^j^^, '^N,j-i,k, ^N~j+i,k^ '^N,j,k-i, and ipN'lk+v we use 
the newly calculated values at the boundaries as soon as they are ready and use the 

f In the remainder of this section, the indices i, j, and k will represent inner grid points (i.e., 
l<i,j,k<N). 
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previous step counterparts for the rest. In practice, ip"-'^ corresponds to the last update 
of such boundaries values. Given the nature of the iterative Crank-Nicholson method, 
these updates occur three times in each time step, making our approach more accurate 
than expected. 

This boundary condition for ip is essential in the implementation of the Hamiltonian 
relaxation method: tests performed using Sommerfeld conditions have shown a 
degradation of the quality of the simulation to BSSN levels. The same tests also show 
that the constraint violating noise generated at the boundaries can be reduced even 
further by generating the expression (|7T|) from the equation 'HN-i,j,k = '^N~2,j,ki instead 
of 'HN-i,j,k = 0, which guarantees a smoother transition of the Hamiltonian residual at 
the next-to-the-boundary points (A^ — l,j,k). 

A problem arises at the edges of the grid, where the expression for the boundary 
value iJN,N,k 

_ F^^r^Ar) Ax Ay 

'PN,N,k ~ n t~ \ ' ^ ' 

^ [lxy)N-l,N-l,k 

is not well defined at the initial time step, when 7^^, = 0. We experimented with several 
different alternatives and found that using a simple bilinear extrapolation of the form 

i'N,N,k = i'N-l,N,k + i'N,N-l,k — i'i,j,k , (23) 

when "jxy falls below some threshold, performs well for more than an orbital period. 



2.3.2. Boundary Conditions for the Gauge Fields a and (3^ The BSSN runs were 
performed using Robin-hke boundary conditions for the gauge fields, to be consistent 
with previous work |2SII2H|- These conditions are based on the asymptotic behavior of 
the gauge fields at distances far from the matter sources jl0| l25]: 

a — 1 oc r^-^ , 

(X y — y u , 
P"^ (XX + X uj , 

(X X y z r^'', 

where the terms proportional to the orbital angular velocity u arise from working in the 
corotating frame {u = u z). 

For the Hamiltonian relaxation runs we found that freezing the lapse function at 
the boundaries to its initial value results in a clear reduction of the incoming constraint 
violating noise. Hamiltonian relaxation test runs were performed using both the /9-freeze 
condition and the F-Driver algorithm with Robin-like boundary conditions. 



2. 3. 3. Hydrodynamical Formulation The evolution of the gravitational fields is coupled 
with the corresponding hydrodynamical equations for the matter sources. For this 
reason, BNS testing is also a valid numerical setup for hydrodynamical formulations and 
their corresponding numerical implementations. In the same way the inspiral phase of 
the BNS offers ideal conditions for the testing of gravitational evolution algorithms, the 
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merger phase can be used to discriminate between different hydrodynamical methods, 
with special interest in their abihty to conserve angular momentum and capture shocks. 
For this paper, we adopt the hydrodynamical formulation and numerical techniques 
that have been used successfully in the simulation of the inspiral phase of BNS I28j. 
since we are only interested in testing gravitational evolution schemes. We describe 
the fluid inside the stars using a perfect fluid stress-energy tensor and a polytropic 
equation of state, with constant F = 2. Section 13.2.21 provides some characteristics of 
the hydrodynamical part of our code. For more details, we refer the reader to the above 
mentioned papers. 

3. Numerical Testing Ground 

We design a testing ground for numerical algorithms that captures some of the main 
characteristics of compact-object binary scenarios: three-dimensionality, full general 
relativistic gravitation, and hydrodynamical evolution. Two types of numerical codes 
are necessary for this endeavor. One to generate astrophysically realistic BNS ID sets 
(CFC-Solver) and another to evolve this data in time (GRHyd). Our codes work on 
Cartesian grids and use finite-difference second order operators within grids that have 
uniform and identical spacing along each axis. The BNS considered in this paper are 
composed of identical stars with equatorial symmetry, to allow for the implementation 
of TT and equatorial symmetry in our codes. We work in the reference frame that rotates 
with the binary, and the stars are aligned with the y axis. 

3.1. Initial Data Code: CFC-Solver 

The core of CFC-Solver is an elliptic solver based on a variation of a multigrid algorithm 
used for previous work [HI |32] • It generates ID sets for BNS through the solution of 
the Hamiltonian and momentum constraints for a quasi-equilibrium circular orbit, via 
the Wilson-Mathews Conformally Flat Condition (CFC) approach mmj. The CFC 
method simplifies the constraint equations by restricting the spatial metric tensor to a 
conformally fiat form and favors "orbit circularity" by imposing a helical Killing vector 
to the spacetime. CFC-Solver can generate data for stars with arbitrary masses and 
spins. For simplicity, the testing ground explored in this article includes only corotating 
and irrotational BNS. We refer the reader to [33] for more details about the code. 

For this article, we use two different initial data sets: one describes a corotating 
BNS in a low-resolution small grid and the other an irrotational BNS in a high-resolution 
large grid. The latter is the same ID set used in previous work |2H], to facilitate the 
comparison of results obtained with our new code. Details of the ID sets are given in 
tabled 
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3.2. Evolution Code: GRHyd 

The results presented in this paper were obtained with a new numerical code for 
the time evolution of gravitational and hydrodynamical fields. GRHyd was written 
using the Cactus programming environment [46 . One of the defining characteristics 
of the code is the separation of the formulations in different and interchangeable 
modules ("thorns" in Cactus language). Each module implements a given set of 
evolution equations, finite-difference schemes, gauge fields, and the corresponding 
boundary conditions. This applies to both the gravitational and hydrodynamical 
formulations. Comparisons between different gravitational formulations are performed 
employing the same hydrodynamical scheme. In this paper we will compare the 
results of the Hamiltonian relaxation scheme with those of BSSN, using in both cases 
a hydrodynamical module based on a van Leer advection scheme jl^ with artificial 
viscosity. 



3.2.1. BSSN and Hamiltonian Relaxation Modules The BSSN module, which follows 
closely the numerical techniques and parameters implemented by Duez et al. jSH] and 
Marronetti et al. ||28j for BNS simulations, is based on a second order in time iterative 
Crank-Nicholson (ICN) scheme. This module plays two essential roles: code validation 
and benchmarking of new formulations. To validate the new code GRHyd, we use the 
BSSN and van Leer modules to perform a simulation based on one of the ID sets used 
in [2H1 (see [Appendix A for the comparison of results). 

Any new numerical method for the evolution of gravitational fields will be matched 
against BSSN simulations. Thus, the BSSN module provides the reference frame on 
which the new schemes will be measured. The first level of comparison involves short 
trial- and-error runs, where poorly performing schemes can be quickly weeded out. After 
this first round, the best performing methods will be used in high-resolution large-grid 
runs, to test their effectiveness in more realistic and demanding simulations. 

The Hamiltonian relaxation module implements the algorithm described in section 
121 The values of the relaxation parameters used in this paper eu = 0.0001 and riu = 70.0, 
were determined empirically. The relaxation is performed until the L2 norm of the 
Hamiltonian residual is smaller than the one at the previous time step for up to a 
maximum of 25 iterations. 

Both BSSN and Hamiltonian relaxation modules use ICN with a Courant factor 
of 0.46. They use the same boundary conditions for the fields 'jij, Aij, K (Sommer- 
feld), and (F* = 0). They also share the same K-Driver and F-Driver parameters: 
e„ = 0.125 t, 7]^ = 0.1, 6/3 = 0.0005, 7]^ = 0.2. The K-Driver (F-Driver) relaxation is 
iterated 5 (10) times. The modules differ, however, in the boundary conditions for the 
lapse: the BSSN module employs Robin-like conditions, while the Hamiltonian relax- 
ation freezes the lapse at the boundaries to its initial value. 



f The value = 0.625 reported in is incorrect. 
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3.2.2. Van Leer Module The van Leer module implements the evolution of the 
hydrodynamical fields using an ICN scheme that couples to the evolution of the 
gravitational fields. The fiuid advection follows a van Leer algorithm [UI that is 
complemented with the introduction of artificial viscosity, to handle the presence of 
shocks. This method, while not as sophisticated as high-resolution shock-handling 
techniques [IH], is fast, easy to program, and good enough for the simulation of the BNS 
inspiral phase as long as no shocks are present. During the inspiral, shocks can develop 
only in the atmospheric envelope that surrounds the stars, and the numerical handling 
of such atmosphere usually requires a good deal of work 39 . We avoid shocks and the 
problems associated with them by adopting the non- atmospheric method introduced in 
We use Copy outer boundary conditions for the hydrodynamical fields. 

3.3. Benchmark Runs 

For the short trial-and-error runs to be practical, the number of grid points has to be 
small enough that the runs can be performed quickly. However, a minimum grid size is 
needed to guarantee simulations of acceptable quality. The short runs are based on a 
corotating (i.e., tidally-locked) BNS system. For the long runs, we use the BNS labeled 
case B in table II of jlHj, which corresponds to an irrotational system outside the ISCO. 
The details for both grids are presented in tabled The number of points in the short run 
grid is about 73 times smaller than that of the long run simulations, making it possible 
to simulate a BNS orbital period in a couple of hours on a typical single-processor 
workstation. The long runs were performed using the IBM p690 Regatta cluster at 
NCSA. 

4. Results 

4.1. Short Trial- and- Error Runs 

In this section we compare BSSN and Hamiltonian relaxation results using short trial- 
and-error runs (see table Q). We concentrate here on the study of the Hamiltonian 
relaxation technique. A more comprehensive study of gauge fields choices and 
boundaries conditions will be done in the future. The BSSN benchmark run uses the K- 
Driver condition for the lapse and the F-Driver for the shift. The Hamiltonian relaxation 
method was tested using the same choice of gauge fields, plus the /3-freeze condition. 
Apart from the use of the Hamiltonian relaxation for determining the conformal factor 
(with its respective boundary condition), the most important difference between this 
method and BSSN is the use of frozen boundary conditions for the lapse as explained in 
sectionEini The results obtained with Hamiltonian relaxation in combination with the (3- 
freeze condition are marginally better than those obtained with the F-Driver algorithm. 
Previous BNS simulations |23 I2H] show that the F-Driver shift condition outperforms 
the /3-freeze condition. In this article, we compare the combination BSSN/F-Driver vs. 
Hamiltonian Relaxation//9-freeze, given that those combinations are the best performers 
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Table 1. BNS Testing: Physical and numerical characteristics of the ID sets and 
Cartesian grids used for testing. The spatial resolution is given in number of grid points 
across the stellar diameter. The bounding box length B gives the extent of the physical 
space covered in each direction (i.e.; the numerical grid spans from[ — B, 0.0, 0.0] 
to [ i?, B, B] since we make use of the equatorial and tt symmetries of the systems), 
in units of total rest mass Mto- The stars where modeled using a polytropic EOS 
with index n=l (F = 2). For this particular EOS, the critical rest mass of a star in 
isolation is mto = 0.180. The binaries are composed of identical stars with individual 
rest masses with 80% of the critical value. The last row gives the number of side-to-side 
light crossing times (let) per orbital period. 





Short Runs 


Long Runs 


Binary type 


corotating 


irrotational 


Total Rest Mass 


0.2920 


0.2932 


Total Grav. Mass Mq 


0.2697 


0.2706 


Orbital ang. vel. fl Mto 


0.03456 


0.02636 


Total Ang. Mom. Jo/M§ 


1.1469 


1.0117 


Spatial Resolution 


20 


40 


Grid Bounding Box B/Afbo 


8.7 


18.6 


Grid points 


3l2 X 60 


128^ X 256 


let per Orbital Period 


10.6 


6.5 



of each formulation. 

Figure Q presents the evolution of the Hamiltonian residual Ti for the first three 
time steps of a BSSN simulation. The same curves are plotted on two different scales, to 
appreciate the effects on the bulk and the boundaries of the grid. We note that, while 
the residual remains unchanged in the bulk of the grid, the violation noise caused by 
Sommerfeld-like boundary conditions grows very fast at the boundaries. This effect is 
accentuated in the short runs, due to the small grid size. Figure El shows the results 
of the Hamiltonian relaxation run. The special boundary conditions of the conformal 
factor ip control very effectively the violation noise generated at the grid's edge. At the 
same time, the Hamiltonian relaxation method reduces the constraint violation in the 
bulk of the grid. 

The effect of the Hamiltonian relaxation technique on the overall quality of the run 
can be better appreciated by looking at the evolution of the total angular momentum 
of the BNS shown in figure El The total angular momentum is the most critical quality 
control curve in binary simulations. Any degradation in the quality of the run is 
appreciated first in the behavior of J. This is the case even for Newtonian simulations 
of binary systems |39j. The angular momentum of the BSSN run is displayed for 
comparison (dashed line). The BSSN run performs reasonably well for about 1/3 of 
an orbital period before it becomes unstable. The Hamiltonian relaxation simulation 
continues, with degrading quality, for over an orbital period. This is a remarkable result 
considering the small size of the short run grid, where an orbital period is roughly 
equivalent to 10 side-to-side light crossing times. 
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Figure 1. BSSN Short Run : Evolution of the HamUtonian constraint residual for 
the first three time steps. The same curves are plotted on two different scales, to 
appreciate the effects on the boundaries (top) and the bulk (bottom) of the grid. The 
curves are plotted following the line with coordinates (0,j/, 0), that runs through the 
center of the star. The companion star, located on the negative y hemisphere, is not 
shown. 




Figure 2. Hamiltonian Relaxation Short Run : Evolution of the Hamiltonian 
constraint residual for the first three time steps. 



Hamiltonian relaxation 



15 



2.00 




Figure 3. Short Runs : Total angular momentum J as a function of time, given as 
fraction of the orbital period P. J is normalized to its initial value Jq. The solid 
(dashed) line corresponds to the Hamiltonian relaxation (BSSN) run. 



4-2. Long High- Resolution Large- Grid Runs 

The Hamiltonian relaxation method was tested with a more accurate simulation, by 
increasing the size and spatial resolution of the numerical grid and starting from an ID 
set corresponding to an irrotational binary (see table Q). We discuss in this section the 
most important features of such simulation, while the corresponding convergence tests 
are covered in [Appendix A All the plots of this section show curves that, for clarity. 



have been normalized to their corresponding initial values. The solid (dashed) lines 
represent the results for the Hamiltonian relaxation (BSSN) run. 

The evolution of the L2 norm of the Hamiltonian constraint residual Ti. across the 
numerical grid is shown in figure IH The satisfaction of the constraint improves by about 
two orders of magnitude with respect to the BSSN run. The suppression of the constraint 
violation is such that, in average, the Hamiltonian residual is five times smaller than the 
ID set value. This drastic reduction of the constraint violation suggests a new way of 
generating ID for binaries: a snapshot of all the gravitational and hydrodynamical fields 
taken aXt/P = 0.5 can be used as an ID set for evolutions with different numerical codes. 
Note that the ID set obtained in this manner would not be subject to the conformal 
flatness restriction for the spatial metric. On the contrary, given that (for this grid) half 
a period is equivalent to more than 3 side-to-side light crossing times, the ID set would 
possess the characteristics (to some extent) of outgoing gravitational radiative fields. 

Figure shows the evolution of the total angular momentum. In addition to the 
Hamiltonian relaxation and BSSN curves, we plotted the PN estimation (dotted line) 
of the angular momentum loss for a point-mass binary with the same mass and angular 
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Figure 4. Long Runs: Evolution of the L2 norm of the Hamiltonian constraint 
violation across the numerical grid. The solid (dashed) line corresponds to the 
Hamiltonian relaxation (BSSN) run. The dotted line marks the violation at the 
initial time step. Note that the curves are plotted in logarithmic scale to highlight the 
more than two orders of magnitude difference between the Hamiltonian relaxation and 
BSSN results. The Hamiltonian relaxation scheme not only suppresses the constraint 
violation modes, but also reduces the violation present in the ID set by a factor of 
about 5. 



momentum as the BNS in consideration (see Appendix B] ). The Hamiltonian relaxation 



run agrees with the PN prediction for about 1.5 orbital periods, while the BSSN curve 
starts going upward well before the first period. The inset of figure El zooms in on the 
first half a period, showing the reduced level of noise in the Hamiltonian relaxation 
curve. 

Figure IHl shows the remaining quality control curves: the coordinate separation 
between stellar centers d, the total gravitational M mass, and the x component of the 
momentum constraint. The total rest mass of the system remains invariant to within a 
0.1 % in both runs. 

We note that the Hamiltonian relaxation run degrades quickly after 3/2 orbital 
periods (about 10 side-to-side light crossing times), while the BSSN run continues for 
over another period before stopping. However, the overall quality of the Hamiltonian 
relaxation results is superior to those of the BSSN simulation during the overlapping 
time. The onset of the instability at the end of the simulation seems to be independent 
on the grid size (see figure IX2I in [Appendix A| ) , possibly indicating that the problem 
is not originated at the boundaries. The total mass and angular momentum of the 
binary are affected much earlier than the Hamiltonian constraint violation, which may 
indicate that the relaxation technique is not at fault either. One possible cause for the 
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Figure 5. Long Runs: Evolution of the total angular momentum J. The solid 
(dashed) line corresponds to the Hamiltonian relaxation (BSSN) run, while the dotted 
line shows the PN estimation (sec [Appendix B| ) . The inset expands the plot for the 
first half orbital period. 
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Figure 6. Long Runs: Remaining quality control curves for Hamiltonian relaxation 
(solid lines) and BSSN (dashed lines). From top to bottom, we show the evolution of 
the coordinate separation between stellar centers d, the total gravitational M mass, 
and the x component of the momentum constraint. 
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premature end of the simulation could be our choice of /3 freezing for the shift vector, 
which is known to perform poorly when the matter distribution drifts significantly from 
the initial configuration. This instability will be more rigorously studied in future work. 

5. Conclusions 

We introduced a new testing ground for numerical relativistic formulations and finite- 
difference techniques, based on short trial-and-error runs of BNS systems. We show the 
usefulness of BNS testing by developing a new technique based of the approximate 
solution of the Hamiltonian constraint (Hamiltonian relaxation). The new method 
was tested using more realistic and computationally demanding numerical simulations 
of BNS outside the ISCO, and the results were compared with those from BSSN 
simulations. The Hamiltonian relaxation improves the overall quality of the simulations 
by reducing the Hamiltonian constraint violation by about two orders of magnitude with 
respect to BSSN values, and by a factor of five with respect to the violation present 
in the initial data set. More remarkably, the improvement in the time evolution of 
the total angular momentum of the binary is such that it agrees with PN point-mass 
binary estimations. The Hamiltonian relaxation run becomes unstable at about 3/2 
orbital periods (~ 10 side-to-side light crossing times), possibly due to the use of (3 
freezing shift condition. This problem will be addressed in future work. Similarly, 
we will attempt to generalize the Hamiltonian relaxation technique to the momentum 
constraint, to further improve the quality of the simulations. Finally, note that the 
Hamiltonian relaxation method can easily be applied to other formulations other than 
BSSN, since only requires the conformal decomposition of the spatial metric tensor. 
This subject will also be explored in the future. 
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Appendix A. Code Tests 

Appendix A. 1. Validation Tests 

Several code validation runs were performed with the new general relativistic code 
GRHyd. Using GRHyd's BSSN and van Leer modules (see section 13. 2p and employing 
the same ID sets, we performed BNS simulations that were compared to those in |28j . 
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Figure Al. Angular Momentum comparison between GRHyd (dashed line) and 
results from Marronetti et al. |2H1 (solid line). The GRHyd simulation starts with the 
long run ID set of table H which corresponds to case B of table II in [2H1 ■ The inset 
zooms in on the first part of the evolution. 

Figure lAll compares the evolution of the total angular momentum of the long run 
simulation described in table ^ with the corresponding case B of table II in |2H] • The 
solid (dashed) line corresponds to the results from (GRHyd). The inset zooms 
in on the first part of the simulation, to highlight the agreement between both codes. 
The differences between the runs are due to two sources: roundoff error, which can 
cause drifts of the order of a fraction of a percent in extensive quantities like the total 
angular momentum, and, more importantly, a grid size difference between both runs. 
The numerical grid employed by GRHyd is one grid zone shorter than the one used in 
[2H]- This is due to the different ways in which both codes setup the symmetries in 
the numerical grid. This effect can be corrected by generating a larger-grid ID set and 
using it for new simulations with both codes. For clarity, we decided to maintain the 
same ID set used in j2H]- The results start diverging significantly from each other after 
1.3 orbital periods. This is expected, since the boundary effects are by then dominant, 
degrading significantly the quality of the simulation. 

Appendix A. 2. Convergence Tests 

We tested the convergence of the results obtained with the Hamiltonian relaxation 
method with the size of the numerical grid. We performed the long run simulation of 
table[T]on four different grids with lengths B/Mho = 9.3, 11.6, 14.0, and 18.6. All grids 
have the same spatial resolution of about 40 points across the stellar radius. The results 
are presented in figure \K2\ where we show the total angular momentum as a function of 
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Figure A2. Convergence of the Hamiltonian relaxation results with varying grid 
sizes. The convergence test is based on the long run simulation (see table ^ and the 
plot shows the behavior of the total angular momentum. The labels for curves are 
(from smallest to largest): dash-dotted, dash-double-dotted, dashed, and solid. All 
the grids have the same spatial resolution. The dotted line shows the PN estimation 
(see [Appendix B| . The inset expands the plot for the first half of the run. 



time. The inset expands the scale for the first half orbital period. The plot shows the 
convergence of the numerical results towards the PN estimation for point-mass binaries. 
The agreement between the three largest grids' results indicates that they could be used 
for preliminary runs, for studies which require extensive parameter space searches, like 
in the determination of the BNS ISCO (see [2H])- 

We also tested the numerical convergence of the Hamiltonian relaxation results with 
the spatial grid resolution. We performed three runs based on the irrotational long run 
of table [U using 20 (low res.), 25 (medium res.), and 40 (high res.) grid points across the 
stellar diameter. Figure IA3I shows the normalized total angular momentum, where the 
separation between curves is consistent with second order convergence in grid spacing. 

In order to see the second order convergence of the Hamiltonian constraint violation 
with grid resolution the relaxation stopping criterion needs to be changed: instead of 
comparing the L2 of the Hamiltonian constraint residual with its value at the previous 
time step (see section EH}, we compare it to the value at the previous iteration step 
(represented by the index m in equation IT^ . The relaxation is stopped when the 
difference between these norms falls below some threshold value. The results for the 
three different grid resolutions are shown in figure IA41 The change in criterion has little 
effect on the results (mass and angular momentum do not vary significantly), however 
the computational time is considerably longer. 
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Figure A3. Total angular momentum in units of total ADM mass squared for the 
irrotational Long Run described in table ^ with three different grid resolutions: low 
(dotted line), medium (dashed line), and high (solid line). The inset shows the ratio 
between the curves (Med — Low) /{High — Med) which, for second order convergence, 
is expected to be about 1.40 (dashed line). 
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Figure A4. Evolution of the L2 norm of the violation of the Hamiltonian constraint 
for the irrotational Long Run described in tablenwith three different grid resolutions. 
Curve labeling is the same as in figure IA3I The numerical factors multiplying the 
curves correspond to the ratios between grid spacings. 
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Appendix B. Post-Newtonian Estimation of Angular Momentum Loss 



In order to get an independent estimation of the angular momentum loss due to the 
gravitational radiation, we calculated the post-Newtonian prediction corresponding to a 
point-mass binary in circular orbit. Kidder [IHj derives the 2.5 PN equations of motion 
for such binaries with arbitrary masses and spins in quasi-circular orbits, which is a good 
approximation when the inspiral motion following the radiation timescale is much larger 
than the orbital motion. For equal-mass zero-spin stars the loss of angular momentum 
becomes 



(Mor) 



1/2 



2423 + 588r] 



336 



Mo 
r 



+ Ait 



Mo 
r 



3/2- 



dJ _ 32 2 [Mo 
dt ~ 5"^ vV 

where Mq is the total mass, rj = M1M2/MI the reduced mass, r the coordinate 
separation, and J the total angular momentum. This formula was derived following 
the Symmetric Trace- Free (STF) radiative multipoles treatment given by Thorne |50jj. 
It requires a formula for the point-mass coordinate separation r as a function of time, 
which can be calculated from the equations of motion as 



dr 
dt 



64 



7] 
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The dotted curves presented in Figs. l5l and lA2l were obtained by integrating these ODEs, 
assuming that the inspiral starts from a stationary circular orbit, with initial values for 
Mq, Jo, and the angular velocity given in tabled The initial coordinate separation 
To is determined by Kepler's law 
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